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Abstract. We investigate a phase-field-crystal model for homogeneous nucleation. 
Instead of solving the time evolution of a density field towards equilibrium we use a 
String Method to identify saddle points in phase space. The saddle points allow to 
obtain the nucleation barrier and the critical nucleus. The advantage of using the 
phase-field-crystal model for this task is its ability to resolve atomistic effects. The 
obtained results indicate different properties of the critical nucleus compared with bulk 
crystals and show a detailed description of the nucleation process. 



1. Introduction 

If a liquid is cooled below its melting temperature the liquid will exist in a metastable 
state until a nucleation event occurs. If the source of nucleation in the undercooled 
melt is only due to fluctuation phenomena the nucleation is called homogeneous. In 
classical nucleation theory a spherical shape for the critical nuclei is assumed and its 
size is determined as a result of competition between the bulk free energy reduction and 
interfacial energy increase. If V is the volume, A the surface area, Ag the chemical 
free energy change per unit volume and 7 the specific interfacial energy, the free energy 
change according to the formation of a new phase is given by AG = V Ag + Ay. For a 
spherical shape with radius r we thus obtain AG = 4/37rr 3 Ag-|-47rr 2 7. The radius r* of 
the critical nucleus must then be such that r* = — 2^/Ag with the critical free energy of 
formation of a critical nucleus AG* = 167r7 3 /(3(Ag) 2 ). This classical theory has been 
utilized to interprete kinetics of many phase transformations and have had some success 
for providing good descriptions on the nucleation kinetics for various systems. On the 
other hand, nucleation is generally significantly more complicated. The shape might not 
be spherical due to an anisotropy of the interfacial energy between a nucleus and the bulk 
phase, which results from the crystallographic nature of a solid nuclei. Furthermore, 
the bulk properties of small nuclei may differ from bulk values typically obtained from 
larger samples. To account for these phenomena various new attempts in the context 
of diffuse interface models have been made to describe nucleation. Such a non-classical 
theory was pioneered by Cahn and Hilliard pQ. For subsequent studies, generalizations 
and specific applications to nucleation, we refer to [21 El IH E] and the references therein. 
In these studies an order parameter is used to distinguish between the nucleus and the 
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bulk phase. Since nucleation takes place by overcoming the minimum energy barrier, 
a critical nucleus is defined as the order parameter fluctuation which has the minimum 
free energy increase among all fluctuations which lead to nucleation. Therefore, the 
critical nucleus can be found by computing the saddle points of the energy functional of 
the order parameter, that has the highest energy in the minimum energy path (MEP), 
which is the path whose highest energy is the lowest among all possible paths. This 
is consistent with the large deviation theory which states that the most probable path 
passes through the saddle point in the large time limit. An efficient numerical approach 
for finding minimum energy paths and saddle points, the so-called string methods (SM), 
has been introduced in [BJ. The method is related to the nudged elastic band (NEB) 
method [7|. Other approach are e.g. the minimax method, which as been used in [8J 
or a phase-field type approach, as used in [U] in the context of nucleation. We will here 
apply a simplified string method (SSM) (TO] but not on an underlying diffuse-interface 
model but a more detailed phase-field-crystal model [11], which accounts for the discrete 
effects on the small length scales involved in nucleation. 

The outline of the paper is as follows: In Sec. 2 we introduce the phase-field-crystal 
model as a local approximation of a classical dynamic density functional theory. In Sec. 
3 we describe the used string method. In Sec. 4 we discuss implementational issues. In 
Sec. 5 we show results for homogeneous nucleation. Conclusions are drawn in Sec. 6. 

2. Phase-field-crystal model 

The phase-field-crystal (PFC) model is by now widely used in order to describe 
solid-state phenomena on atomic length scales. The PFC model was first developed 
in [TT] and then subsequently applied to many situations like interfaces [121 [IB"] , 
polycrystalline pattern formation [T4"l [15], commensurate-incommensurate transitions 
[IB] , edge dislocations [17], grain boundary pre-melting [IS] , colloidal solidification 
[19] and dislocation dynamics [20J. The model resolves the atomic-scale density wave 
structure of a polycrystalline material and describes the defect-mediated evolution of 
this structure on time scales orders of magnitude larger than molecular dynamic (MD) 
simulations. In its simplest form the PFC model results from the energy 

F[<p]= J^-\W V \ 2 + \^ V f + f^)dx (1) 

with /(</?) = |(1 — e)f 2 + |y? 4 a potential, ip the number density and e a parameter 
determining the approximation of the liquid structure factor [II] . Comparing the 
energy with a classical phase-field type energy, e.g. J n ||V0| 2 + \g(4>) dx for an order 
parameter 0, with 5 a length scale determining the width of a diffuse interface and g(4>) a 
double well potential, the difference is in the sign of the gradient term and the additional 
higher order term. The negative sign in the gradient term favors a changes in ip, whereas 
the higher order term favors to suppress such changes. The competition between both 
terms introduces a fixed length scale for which the energy will be minimized. This length 
scale is used to model the periodicity of a crystal lattice. The formulation used here 
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favors a hexagonal closed packed structure in two dimensions. Due to the underlying 
periodicity, several solid state phenomena such as elasticity, plasticity, anisotropy and 
multiple grain orientations are naturally present in the formulation. The dynamic law 
constructed to minimize the free energy follows as the if ~ ^gradient flow of the energy 



dt<p = A/i (2) 

Sip 



with chemical potential pi = and the variational derivative given by 



dip 

This defines the PFC model and its evolution is by construction towards a (meta) stable 
state. 

Although this formulation is phenomenologically, the PFC model can be derived 
starting from a Smoluchowski equation via dynamic density functional theory using 
various approximations [2TJ HH] • With an appropriate parameterization it thus provides 
also a quantitative atomic theory, operating on diffusive time scales. Within new 
developments [TH [23] quantitative agreement in computed properties could be achieved 
using the PFC model for various materials. It thus provides an ideal model to study 
nucleation. 



3. Minimum energy path 

3.1. Definition 

The dynamics shown in eq. ^ describe the evolution towards equilibrium of a single 
state tp in phase space according to the generalized thermodynamic force A/i. But in 
order to characterize nucleation the most likely transition path between (meta) stable 
states has to be identified. In the description of such a non equilibrium process, the 
minimum energy path (MEP) plays a crucial role. The MEP is a path in phase space 
that connects (meta) stable states. A path in phase space is thereby defined as 

7 C = {y?« : a e [0,1]} 

with a a parameterization of the path. For the MEP the generalized thermodynamic 
force Afi is tangential to this path: 

(A//) 1 = (3) 

Thus, using the dynamics in eq. ([2]), any state of the MEP will always evolve along the 
MEP towards a stable state. That is, the MEP is a real reaction path in phase space 
and along the MEP the energy is defined by eq. ([!]). The local energetic maxima and 
minima along the MEP can be used to determine the nucleation barrier and the critical 
nucleus. The string method (SM) is been designed to find the MEP. It evolves a given 
chain of states towards the MEP, see Fig. [T] for illustration. 
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Figure 1. Schematic sketch of a free energy surface of a two dimensional 
phase space. Two (meta) stable states (valleys) are separated by a saddle 
point. The single states (circles) are evolved by thermodynamic forces 
(arrows). The initial string (straight line) is evolved by the string method 
towards the MEP. 



3.2. String method 

A path in phase space 7 C is represented by a discrete set of states (fi, denoted by 

T = { Vi ■ i = o, 1, .. .,N} = {(ft} . 
The length of the path is defined by 



N-l 



L d) = - and m = / i^ r )i dr 

i=0 J 

where | • | measures the distance between two states and is defined by the L 2 -norm. Thus, 
the normalized tangent i to the path at state (fi may be calculated as ii = 

The idea behind the SM is to minimize the free energy of all single states according 

to 

ipp 1 = rA/i + . (4) 

but restricting the evolution to orthonormal direction of the path. In addition a 
tangential force is included in order to keep the quality of the representation of the 
path by the string. 

7 n ^7™ +1 =W- +l ) (5) 
with < +1 = r (A/i)^ + < + \t 

The Langrang multipliers \ are e.g. uniquely determined by forcing an equidistant 
distribution of states along the path, r is a fictitious timestep and controls the velocity 
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of evolution. Eq. ^ defines the String Method. It is easily seen that the MEP 
7mep is an invariant according to the dynamics of SM. By definition of the MEP, the 
thermodynamic force is only tangential to the path. The constraints introduced by the 
tangential force do not alter the path, but only reparameterize the representation of 7. 
Thus, 7mep an d 7mep represents the same path in phase space. 

In order to implement the SM the thermodynamic force has to be calculated 
and projected to the orthogonal direction of the path. Furthermore the Langrange 
multipliers have to be calculated. In order to simplify the calculation the SM can be 
divided in two steps leading to a Simplified String Method (SSM). First the string is 
evolved due to the thermodynamic force and then the path is reparameterized, see [TD] . 
That is, the states representing the path {(pi) are replaced by equally distant states, 
that represent the same path {(fi\. This new states are constructed by interpolation 
between the original states {<fi}. As in every evolution step there is a parameterization 
step, it is no longer necessary to project the thermodynamic force. 

Thus the SSM is defined by two steps: 

(i) Evolution step: 

i n ^r +1 = {tf +l ) (6) 

with = tA V + <Pi 

(ii) Reparameterization step: 

7 n+1 ^7 n+1 = K +1 } (7) 



suchthat|< + + 1 1 -^ +1 | = ^tj f * = 0,l,2,...,iV-l 
and 7 ra+1 and 7 n+1 representing the same path in phase space. 

Here the reparameterization is done to force equidistant states on the path. However, 
the reparameterization may also be changed to account for problem specific details, e.g. 
to get finer representation at the saddle point. The advantage of SSM over SM is that 
the thermodynamic force has not to be projected. Additionally it is shown that this 
modification leads to a more stable and accurate method |10|. 



L(7 



3.3. Fixed length Simplified String Method 

In the above defined SSM an initial path in phase space is evolving towards the MEP. 
The first and the last state representing the path thereby converge to different (meta) 
stable states. The saddle point or here the critical nucleus is defined by the state of 
highest energy in the MEP. If there is only one energetic maximum and the saddle point 
is well defined, the MEP could be calculated easily by just solving the time evolution 
of a small perturbation of the critical nucleus towards the stable states according to 
eq. (]2]). Thus, only two time dependent simulations have to be done. Therefore, it is 
enough to find the saddle point to reconstruct the whole MEP efficiently In order to 
concentrate the simulation effort to find only the saddle point, we introduce a Fixed 
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Length Simplified String Method (FLSSM). The total length of the string is restricted 
by the reparameterization step. Assume that we allow a maximal length of the string, 
-kfixed- For simplicity we also assume that the first state converges towards a stable state. 
Then, the reparameterization can always project the states back to a string beginning 
with the first state with length Lfi xe d- The last state is not converging to a meta stable 
state, but might be some unstable state but within a different basin. This state can be 
used to reconstruct the whole MEP by a time evolution according to eq. ([2]). The same 
idea can be used at the same time on both sides of the saddle point, by fixing a state 
near the saddle point and restricting the length on both sides of the path. As we need 
some information about the string length and the position of the saddle point, we use 
the FLSSM in order to refine and improve accuracy of a MEP which was calculated by 
the standard SSM with only a few states. The method can be viewed as an adaptive 
approach which efficiently finds the saddle point within a given tolerance. 



3.4- Implementation of a Fixed Length Simplified String Method for PFC 

The string is a set of density distribution {<Pi}. As we consider a closed system and have 
mass conserving dynamics, we have to restrict the possible states representing a string 
to the same mean density, (p — J (pi(r) dr for all i. 

For every state the standard dynamics has to be solved according to eq. §2§. The 
partial differential equation of 6th order is splitted into a set of three second order 
equations: 

d t <pi = Afi 

fi = Av + 2A(fi + f'((pi) 
v = Aipi 

for which a stable semi-implicit finite element discretization is introduced in [22]. We 
use this approach but with higher order elements. The algorithm is implemented in the 
adaptive finite element toolbox AMDiS [21]. 

The fictitious timestep is adjusted such that the reparameterization step can be 
done mostly considering only neighbouring states and such that the evolved state is 
substantial different from the previous one. 

In this work we use linear interpolation between the states. We define the length 
of the string up to state M in analogy to L{p() for the whole string by La/ (7) = 
YI^Lq 1 Wi+i ~ ViY The distance between states after reparameterization is I = 
Then the reparameterized state (pi at L* = il is constructed by linear interpolation 
using the neighboring states of (pi from 7. 

fi = <Pk + (<fk+i ~ fk)ct (8) 

a= fl ~ Lfe( _ 7) . and L k tf) < L* < L k+1 (j) (9) 
\fk+i - <Pk\ 

The FLSSM is implemented in a parallel way. That is, every state define a process 
and the evolution step is calculated in parallel. The result is then send to the nodes 
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of the neighbouring states. In order to avoid complicated communication between the 
processes, the reparameterization step is further simplified. The linear interpolation 
described above is used if the reparameterized state is in between the neighbouring 
states. If not, the reparameterized state is set to one of the neighbouring sites. 



This does not alter the MEP but only the dynamic of the string in phase space, as long 
as we choose the fictitious timestep in a way that for pathes near the MEP at the end 
of simulation always linear interpolation between neighbouring sites can be used. 

We define convergence of the method if the string changes only in tangential 
direction. This is ensured by two criteria. First, the change in string length and the 
change in energy in every state after reparameterization is below a given tolerance which 
ensures no change of the string in normal direction. The second criteria accounts for 
the change in tangential direction and ensures the evolution towards the (meta) stable 
states. Therefore the fictitious time step is adjusted to evolves the states substantially 
far but still allowing for linear interpolation with 0.3 < a < 0.7. If the second criteria 
is not fullfilled, the fictitious time step is adapted. Due to small thermodynamic forces 
near the saddle point, the second criteria may be relaxed for states very close to the 
saddle point. 

4. Results 

We consider as a proof of concept the nucleation of a crystal grain in an undercooled 
liquid. The parameters needed in the PFC model, eq. (|2]), are the mean value of the 
density, (p — J (p(r) dr and the parameter e, which can be interpreted as a driving 
force of the phase change, e.g. undercooling [261 125] or strength of interaction [19]. 
In our example we choose parameters in the coexistence region of the phase diagram 
(e, (p) = (—0.289,-0.345). For this parameter the liquid is a meta stable state. The 
stable state is a grain in coexistence with liquid. The grain is slightly anisotropic and 
there is a small density difference between crystal and liquid, see [13] . 

In order to calculated the MEP an initial string 7 = {$} has to be defined 
such that the mean density of every state is equal <p> — (p® and that the first and the 
last state evolve towards two different (meta) stable states. In our work, we use two 
different initial strings to demonstrate that the obtained MEP is independent of the 
initial configuration. In the first example we set the first state to liquid ipo(x) = (p 
and the last state to the equilibrium shape of the grain. The states in between are 
then constructed by linear interpolation. In the second example, every state was set 
homogeneous and disturbed by white noise r\ such that <Pi(x) = (p + SiT). The strength of 
the noise Si, was linearly increased starting from to represent the liquid state towards 
a large value, which ensures evolution within the coexistence regime. Both initial strings 




(10) 
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converge to the same MEP shown in Fig. [2] In order to proof stability of the obtained 
MEP every state was disturbed independently by some random field and than taken as 
a initial string to recalculate the MEP. 




-3.5 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 

normalized length of string 



Figure 2. Free energy along the MEP. The free energy is plotted relative 
to the liquid state at normalized string length, 1=0. The path is discretized 
by 94 states. At 1=1 the stable state is reached. 

The first state corresponds to liquid and the last to a grain in coexistence with the 
liquid. The grain equilibrium is energetically favourable compared to the liquid and is 
the stable state in phase space. The liquid state is meta stable. The nucleation barrier 
or the saddle point is found at normalized string length of approx. 0.04. States right 
to the saddle point correspond to growing crystallites and left to melting crystallites. 
The string was discretized by 94 states which are equally distributed, so the region 
around the nucleation barrier is resolved only by 10 states. In order to get a better 
resolution the FLSSM is used. The length of the string is therefore restricted to | of the 
original length of the MEP and is rediscretized by 46 states, which are constructed by 
linear interpolation of the calculated MEP. We can view this as an adaptive method to 
increase the accuracy of the calculated saddle point or a proof that the obtained saddle 
point is independent of the used parameterization of the string. In our example this 
independency is shown. Fig. [3] shows the obtained nucleation barrier AE and critical 
nucleus. 

The critical nucleus is defined by the state indicated by (c), <p c . (b) indicates y?b 
a sub critical state, which most likely will melt, (d) - (f) indicate states (fd - (fif which 
correspond to super critical states which will solidify. 

In Fig. [4] the density field of the labeled states are shown. The critical nucleus 
is a hexagonal cluster with only seven maxima. A small perturbation of this state 
will lead either to growth towards the equilibrium shape or to melting. The grow is 
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Figure 3. Detailed MEP around the saddle-point. Closed symbols indicate 
the MEP calculated by SSM as in Fig. [2j Open symbols show the result 
achieved by restricting the string length to | and using FLSSM. 



symmetric and can be seen more quantitative in Fig. [5j which shows the density profile 
along the x-axis in the various states of the growth process. The density plot shows 
that the maximum amplitude of the critical nucleus is smaller than in the final bulk 
state. This can correspond to defects in the crystal, as we consider here only a mean- 
field description, or weaker ordering of particles. In both cases this shows that the 
critical nucleus has different structure and bulk energy than the corresponding bulk 
state. Nucleation thus begins with a disturbance that reflects the crystal structure but 
has a small amplitude. During growth the spatial size of the initial fluctuation and the 
amplitude increases. At state (p e the maximum amplitude is equal to the bulk value 
and does not increase anymore. After this state the grain begins to grow only along the 
solid liquid phase boundary. 

These computations indicate that the subcritical and supercritical grains around the 
critical nucleus are different to ideal bulk crystals. In addition to the non-spherical shape 
of the nucleus, which is not considered in classical nucleation theory, such size dependent 
properties are also not considered in full detail in classical phase field approaches for 
nucleation. 

5. Conclusion 

A phase-field-crystal model is used to determine nucleation barriers and the critical 
nucleus in homogeneous nucleation. The results obtained indicate details of the 
nucleation process which are not considered in classical nucleation theory but also cannot 



Critical nuclei 



10 





Figure 4. Various states at the MEP. The labels (b)-(f) indicate the 
different state tpt, - tp?. 




2 4 6 

length in lattice const 



Figure 5. Density profile of grains at the MEP. The labels (b)-(f) indicate 
the different state <p a - (pt, see Fig. [3] 
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be addressed in full detail with classical phase field models. The obtained size of the 
critical nucleus, which here consists only of 7 atoms furthermore asks for an atomistic 
description. Even if only phenomenological values are used in the computation, the 
described method gives a proof of concept. The described String Method is independent 
of the parameterization of the underlying evolution model and thus will allow also to be 
used for specific materials. With the described implementational details of the String 
Method, concerning parallel processing and adaptive concepts we believe the approach 
to be applicable also in three dimensions. 
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